package investigation;

import java.util.Random;

import precision.Arbitrary;
import precision.Doubles;
import precision.FPObj;
import precision.Floats;
import precision.Ints;

import mathtools.Matrix;

public class PowerMethod {

	/**
	 * Power iteration algorithm
	 * @param a Matrix A
	 * @param e Initial guess for eigenvector e (none zero elements please)
	 * @return Guess of the eigenvalue
	 */
	public static FPObj powerIteration(Matrix a, Matrix e, int percent)
	{

		FPObj ev = (a.multiplyMatrix(e)).getElement(0, 0).divide(e.getElement(0, 0));
		FPObj evold;
		FPObj dev;
		int count = 0;
		int maxcount = 200;
		do
		{
			evold = ev;
			e = a.multiplyMatrix(e);
			e.divideMatrix(e.columnNorm());

			ev = (a.multiplyMatrix(e)).getElement(0, 0).divide(e.getElement(0, 0));
			dev = ev.sub(evold);
			if ( count > maxcount)
			{
				System.out.println("Break");
				break;
			}
			count++;
		}
		while ((dev.abs()).divide(ev).greaterThan(ev.newFPObj(1).divide(ev.newFPObj(100*percent))));
		return ev;
	}

	/**
	 * @param args
	 */
	public static void main(String[] args) {

		// Square size of matrix
		int N_max = 500;


		for(int N=100; N<N_max;N+=50){

			System.out.println("N: " + N);

			int counter = 0;

			for(int B=0; B<50;B++){

				// Random class
				Random rand = new Random();

				// Sparseness of our matrix (i.e. 0-1 chance of getting zero element
				double sparseness = 0;

				// Percent we wish the solution not to change between iterations (i.e. 0.001 %)
				double percent = 0.00001;

				// Array of our precision objetcs
				FPObj[] precs = {new Floats(), new Doubles(), new Arbitrary(0, 20)};

				// Generate the lower diagonal matrix
				int[][] A = new int[N][N];
				for ( int i = 0; i < N; i++ )
					for ( int j = 0; j < N; j++ )
					{

						if ( rand.nextFloat() > sparseness)
						{
							A[j][i] = rand.nextInt(1000);
						}
						else
						{
							A[j][i] = 0;
						}

					}


				// Convert the matrix to Matrix class
				Matrix Aints = new Matrix(N, N, new Ints());
				for ( int i = 0; i < N; i++ )
				{
					for ( int j = 0; j < N; j++ )
					{
						Aints.setElement(i, j, new Ints(A[i][j]));
					}
				}


				// Convert A and V to our precisions
				Matrix[] Ap = new Matrix[precs.length];
				Matrix[] Vp = new Matrix[precs.length];
				for ( int k = 0; k < Ap.length; k++)
				{
					Ap[k] = new Matrix(N,N,precs[k]);
					for( int i = 0; i < N; i++ )
					{
						for ( int j = 0; j < N; j++ )
						{
							Ap[k].setElement(i, j, precs[k].newFPObj(((Ints)Aints.getElement(i, j)).getValue()));
						}
					}
					Vp[k] = new Matrix( N, 1, precs[k]);
					for( int i = 0; i < N; i++ )
					{
						Vp[k].setElement(i, 0, precs[k].newFPObj(1));
					}
				}


				FPObj[] e = new FPObj[precs.length];

				for ( int i = 0; i < precs.length; i++)
				{

					e[i] = powerIteration(Ap[i],Vp[i], (int)Math.round(1.0/percent) );
					//System.out.println("Eigenvalue using precision: " + precs[i].toType());
					//System.out.println(e.toString());


				}
				


				boolean check = false;
				for(int i = 1; i< precs.length; i++){
					for(int j = 0; j<i; j++){	
						if(i==j){
							//Do nothing
						}
						else{
							if((int)(e[i].toFloat())!=(int)(e[j].toFloat()) && !check){
								check=true;

							}
						}
					}
				}

				if(check){
					counter += 1;
					System.out.println("Count: "+counter);
					for(int i=0;i<precs.length;i++){
						System.out.println(e[i].toFloat() + " " + e[i].toType()); 
					}
				}

			}
		}
	}

}
